import math

GP2 = (math.sqrt(5) - 1)/2
GP1 = 1 - GP2


def solve_max(func, a, b, eps=1e-6):
    """
    求解函数的最大值及其对应的解
    :param func 被求解最大值的多元函数
    :param a 每一元的下界。对N元函数来说，a必须是含有N个元素的元组或列表
    :param b 每一元的上届。对N元函数来说，b必须是含有N个元素的元组或列表
    :param eps 当两个可行解之间的距离小于eps时，求解会中止。
    :return 返回二元组x, value，value是最大值，x是对应的解
    """
    if type(a) not in (list, tuple):
        assert type(b) not in (list, tuple)
        a = [a]
        b = [b]
    else:
        assert len(a) > 0
        assert len(a) == len(b)
    num = len(a)
    x = [0] * num
    try:
        value = _solve(x, 0, func, a, b, num, eps)
    except _BreakException as e:
        x, value = e.x, e.value
    return x, value


def _solve(x, i, func, a, b, num, eps):
    if i >= num:
        return func(*x)

    ai = a[i]
    bi = b[i]

    g1 = ai + GP1 * (bi - ai)
    g2 = ai + GP2 * (bi - ai)
    x[i] = g1
    v1 = _solve(x, i+1, func, a, b, num, eps)
    x[i] = g2
    v2 = _solve(x, i+1, func, a, b, num, eps)
    while abs(ai - bi) > eps:
        if v1 > v2:
            bi = g2
            g2 = g1
            v2 = v1
            g1 = ai + bi - g2
            x[i] = g1
            v1 = _solve(x, i+1, func, a, b, num, eps)
        else:
            ai = g1
            g1 = g2
            v1 = v2
            g2 = ai + bi - g1
            x[i] = g2
            v2 = _solve(x, i+1, func, a, b, num, eps)
    return (v1 + v2) / 2


class _BreakException(Exception):
    def __init__(self, x, value):
        super(_BreakException, self).__init__()
        self.x = x
        self.value = value


def break_solve(x, value):
    raise _BreakException(x, value)


if __name__ == '__main__':
    import time
    def f(x1, x2, x3):
        return 3 - (x1 - 2)**2 - (x2 + 3)**2 - (x3/2 -5)**2
    start = time.time()
    x, value = solve_max(f, [-100, -100, -100], [100, 100, 100])
    print('======= test1: =======')
    print('x =', x, '\nvalue=', value, '\nmilliseconds expired:', (time.time() - start) * 1000)

    def f(x1, x2, x3):
        value = 3 - (x1 - 2)**2 - (x2 + 3)**2 - (x3/2 -5)**2
        if value > 1.0:
            break_solve((x1, x2, x3), value)
        else:
            return value
    start = time.time()
    x, value = solve_max(f, [-100, -100, -100], [100, 100, 100])
    print('\n======= test2: =======')
    print('x =', x, '\nvalue=', value, '\nmilliseconds expired:', (time.time() - start) * 1000)
